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The eigenmode spectrum is a fundamental starting point for the analysis of plasma stability and 
the onset of turbulence, but the characterization of the spectrum even for the simplest plasma 
model, ideal magnetohydrodynamics (MHD), is not fully understood. This is especially true in 
configurations with no continuous geometric symmetry, such as a real tokamak when the discrete 
nature of the external magnetic field coils is taken into account, or the alternative fusion concept, 
the stellarator, where axisymmetry is deliberately broken to provide a nonzero winding number 
(rotational transform) on each invariant torus of the magnetic field line dynamics (assumed for 
present purposes to be an integrable Hamiltonian system). 

Quantum (wave) chaos theory provides tools for characterizing the spectrum statistically, from 
the regular spectrum of the separable case (integrable semiclassical dynamics) to that where the 
semiclassical ray dynamics is so chaotic that no simple classification of the individual eigenvalues is 
possible (quantum chaos). 

The MHD spectrum exhibits certain nongeneric properties, which we show, using a toy model, 
to be understable from the number-theoretic properties of the asymptotic spectrum in the limit of 
large toroidal and poloidal mode (quantum) numbers when only a single radial mode number is 
retained. 

Much more realistically, using the ideal MHD code CAS3D, we have constructed a data set of 
several hundred growth-rate eigenvalues for an interchange-unstable three-dimensional stellarator 
equilibrium with a rather flat, nonmonotonic rotational transform profile. Statistical analysis of 
eigenvalue spacings shows evidence of generic quantum chaos, which we attribute to the mixing 
effect of having a large number of radial mode numbers. 

PACS numbers: 52.35.Bj,05.45.Mt 
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I. INTRODUCTION 



The tokamak and stellarator fusion concepts both seek 
to contain a plasma in a toroidal magnetic field, but to 
a good approximation the tokamak field is axisymmet- 
ric. Thus the study of the spectrum of normal modes 
of small oscillations about equilibrium is simplified by 
the existence of an ignorable coordinate. The stellarator 
class of device, on the other hand, is inherently nonax- 
isymmetric and the lack of a continuous symmetry means 
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there are no "good quantum numbers" to characterise the 
spectrum. 

This makes the numerical computation of the spectrum 
a challenging task, but numerical matrix eigenvalue pro- 
grams, such as the three-dimensional TERPSICHORE 
lU and CAS3D codes, are routinely used to assess 
the ideal magnetohydrodynamic (MHD) stability of pro- 
posed fusion-relevant experiments. An example is the 
design of the 5-fold-symmetric Wendelstein 7-X (W7-X) 
stellarator, where CAS3D was used |3j to study a number 
of different cases. 

The configurations foreseen to be studied experimen- 
tally in W7-X are MHD stable, but in the present paper 
we are concerned with an unstable case from this study, a 
high-mirror-ratio, high-rotational transform equilibrium 
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FIG. 1: Left: plasma boundary of a 5-periodic toroidal 
equilibrium geometrically related to W7-X configurations (at 
(/3) = 0.05, VMEC calculation). Right three frames: vertical 
cross-sections of the configuration at three different toroidal 
angles within a field period, with Z the vertical coordinate 
and R the distance from the Z axis. The plasma boundary 
has been scaled such that the minor radius of the torus a ~ I 
and the major radius equals the aspect ratio R/a = A. 



(Fig. r2J- Due to its less pronounced shaping this case 
is quite unstable, which contrasts with the properties of 
genuine W7-X configurations. The three-dimensional na- 
ture of the equilibrium breaks all continuous symmetries, 
coupling both poloidal (to) and toroidal (n) Fourier har- 
monics and thus precluding separation of variables and 
simple classification of the eigenvalues. 

These eigenvalues, lo 2 = —-f 2 , are real due to the 
self-adjointness of the linearized force and kinetic en- 
ergy operators in ideal MHD linearized about a static 
equilibrium. This is analogous to the Hermitian nature 
of quantum mechanics, so we might a priori expect to 
be able to take over mathematical techniques used in 
quantum mechanics. Thus we study the W7-X Mercier 
(interchange)-unstable case mentioned above using sta- 
tistical techniques from the theory of quantum chaos |3, 
eg]. 

This is of practical importance for numerical analysis 
of the convergence of eigenvalue codes because, if the sys- 
tem is quantum-chaotic, convergence of individual eigen- 
values cannot be expected and a statistical description 
must be used. However, there is a fundamental question 
as to whether the ideal MHD spectrum lies in the same 
universality class as typical quantum mechanics cases. 

This question has been addressed recently 0, by 
studying the interchange unstable spectrum in an effec- 
tively cylindrical model of a stellarator. In the cylindrical 
case the eigenvalue problem is separable into three one- 
dimensional eigenvalue problems, with radial, poloidal, 
and toroidal (axial) quantum numbers I, to, and n, re- 
spectively. If the spectrum falls within the generic quan- 
tum chaos theory universality class for integrable, non- 
chaotic systems |8( then the probability distribution func- 
tion for the separation of neighboring eigenvalues is a 
Poisson distribution. 

However, this work indicates that the universality class 
depends on the method of rcgularization (ie truncation 
of the countably infinite set of ideal-MHD interchange 
growth rates): a smooth, physically-motivated finite- 
Larmor-radius roll-off in the spectrum appears to give 
the generic Poisson statistics for separable systems, but 
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FIG. 2: Number N(E) of renormalized eigenvalues E m ,n be- 
low a given value E for m max = 10: (a) for the Hamiltonian 
H — pejy/pe; (b) for the Farey sequence. 

a sharp truncation in to and n gives highly non-generic 
statistics. The latter case is less physical, but corre- 
sponds closely to the practice in MHD eigenvalue studies 
of using a restricted to, n basis set but a relatively fine 
mesh in the radial direction. 

A careful analysis of the spectrum of ideal-MHD in- 
terchange modes in a separable cylindrical approxima- 
tion revealed non-generic behaviour of the spectral 
statistics — a bimodal PDF, rather than the expected 
Poisson distribution. The non-genericity of this separa- 
ble case indicates that caution must be applied in apply- 
ing conventional quantum chaos theory in non-separable 
geometries. 

The study indicated that the non-generic behaviour 
of ideal-MHD interchange modes was due to the peculiar 
feature of the dispersion relation for these modes that the 
eigenvalues in the short-wavelength limit depend only on 
the direction of the wave vector, not on its magnitude. 
(This is unusual behaviour, but it is shared with inter- 
nal gravity waves in geophysical fluid dynamics.) It was 
suggested in Q that the detailed features of the spec- 
trum could be understood from the properties of Farey 
sequences. 

In the present paper we discuss fictitious eigenvalues 
from a toy model, used to elucidate the importance of 
number-theoretic effects, that illustrate how nongeneric 
the MHD eigenvalue spectrum can be if only one radial 
eigenmode is used. Then we present the results of the 
quantum chaos analysis of the W7-X case. 

II. TOY EIGENVALUE PROBLEM 

To gain insight into the nongeneric behavior found in 
the separable case with only one radial eigenmode we 
study the energy spectrum {£ n ,m} for a "toy" quantum 
mechanical Hamiltonian H = p^/pe where the configu- 
ration space is the 2-torus 9 € [0,2-7r), <f> £ [0, 27r) with 
periodic boundary conditions. In the semiclassical ap- 
proximation we see that H depends only on the direction 
of p, not its magnitude, as for MHD interchange modes 
and internal gravity waves. 

The eigenvalue problem is the time-independent 
Schrodinger equation, Hip = Sip, the eigenf unctions be- 
ing exp[i(TO(9 + n</>)]/47r 2 , where to and n are integers. 
The eigenvalues £ n ,m = n/m (to ^ 0). 
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Note the singular nature of the spectrum — it is dis- 
crete, yet infinitely dense, the rationals being dense on 
the real line. Also, the spectrum is infinitely degenerate 
as eigenvalues are repeated whenever m and n have a 
common factor. Mathematically such a spectrum, nei- 
ther point nor continuous, belongs to the essential spec- 
trum |9j. 

In order to analyze this spectrum using standard quan- 
tum chaos techniques we first regularize it by bounding 
the region of the to, n lattice studied, and then allowing 
the bound to increase indefinitely. Fortunately the PDF 
P{s) is independent of the precise shape of the bounding 
line when we follow standard practice [f| in renormalizing 
(unfolding) the energy levels to make the average spac- 
ing unity. Thus we adopt the simplest choice, taking the 
bounded region to be the triangle < n < m < TO max . 
As the points (n, to) form a lattice in the plane with mean 
areal density of 1, we can estimate the asymptotic, large- 
TO max behaviour of the number of levels, iV max , from the 
area of the bounding triangle: the to axis, the line n — m 
and the line to = 1, which gives the "Weyl formula" 

^max ~ TO max/2- 

The Farey sequence T(Q) is the set of all rational num- 
bers p/q between and 1, 1 < q < Q, arranged in order 
of increasing arithmetic size and keeping only mutually 
prime pairs of integers p andg. Farey sequences are im- 
portant in number theory [lOj and have application in 
various dynamical systems problems, such as the the- 
ory of mode locking in circle maps . They even have 
a connection with the famous Riemann hypothesis |12| . 
They are conveniently generated using the Mathematica 
algorithm ^| 

Farey [Q_] := Union [{0,1}, 

Flatten [Table [p/q, {q, Q},{p, q-1}]]]. 

The list Q(m max ) = {£n,m}, sorted into a non- 
decreasing sequence {£i\i = 1,2,... , -/V max } is very sim- 
ilar to the Farey sequence J^(TO m a X ) except for the high 
degeneracy (multiplicity) of numerically identical levels, 
especially when n/m is a low-order rational. 

Define the renormalizcd (unfolded) energy as E n:in = 
Nma.x£n,m- The normalization by iV mmax ensures that 
£jf mi> . = N maX) so the mean slope of the Devil's staircase 
shown in Fig.^a) is unity. The large vertical steps visible 
in Fig. |2f a) are due to the high degeneracy at low-order 
rational values of n/m. 

This high degeneracy is also the cause of the delta 
function spike at the origin visible in the level-separation 
probability distribution plot shown in Fig.[3Ja), which is 
very similar to Fig. 9(a) in 0. 



A. Farey statistics 

The tail in Fig.^a) is due to the non-degenerate com- 
ponent of the spectrum, obtained by reducing all frac- 
tions p/q to lowest terms and deleting duplications. Thus 
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FIG. 3: Separation statistics: (a) for the model Hamiltonian; 
(b) for the Farey sequence. The solid curves are from the 
Farey spacing measure, Eq. In (b) the short-dashed curve 
is for the Poisson Process of the generic integrable problem 
and the long-dashed curve is that for the Gaussian orthogonal 
ensemble of random matrices (quantum chaotic case). 



the eigenvalues in this set are Nq times the terms of the 
Farey sequence T{Q). 

To study the statistics of this non-degenerate compo- 
nent it is natural to define the Farey spectrum {Ef} as 
N F (Q) times the terms of the Farey sequence J-(Q), 
where N F (Q) is the number of terms in J-(Q). The 
asymptotic behaviour of N F (Q) in the large- n limit is 
given [3 p. 391] by N F (Q) ~ 3Q 2 /vr 2 + O(QlnQ). 
The staircase plot and separation distribution P F (s) for 
the Farey spectrum are given in Fig. Etb) and Fig.[3Jb), 
respectively. 

It is a standard result in the theory of Farey se- 
quences 0, p. 301] that the smallest and largest nearest- 
neighbour spacings in J-(Q) are given respectively by 
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so that the support of the tail component of P(s) in 
Fig- EH a ) becomes [1/2, oo) in the limit Q — > oo, while 
that of P F (s) in Fig. GJb) is [3/tt 2 , oo). 

Augustin et al. eq. (1.9), derive the spacing den- 
sity for the Farey sequence as 
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The solid curve in Fig. EHb) is obtained by setting 
P F (s) = gi(s) and is seen to agree well with the nu- 
merical results. 

The solid curve in Fig. |3{a) is obtained by setting 
P(s) = [N F (Q)/N(Q)} 2 gi (N F s/N) and agrees well with 
the tail of the histogram. The ratio of the area of the tail 
in Fig.^a) to the strength of the delta function Fig.^a) 
is N F (Q)/[N(Q) - N F (Q)} « 1.55. 

We have verified that the probability distributions re- 
main unchanged if subranges of the spectra are used , in 
agreement with the result included in Theorem 1 . 1 of [l2| 
that the convergence to a probability measure is indepen- 
dent of the interval chosen. 
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FIG. 4: Left: the rotational transform * = *(s) profile versus 
normalized toroidal flux s (s oc r 2 near the magnetic axis). 
Right: a measure of the Mercier stability versus normalized 
toroidal flux for the 3-dimensional W7-X-like case of Fig. 
A negative value indicates instability. 
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FIG. 5: Choice of basis set of toroidal and poloidal Fourier 
harmonics. 



III. W7-X RESULTS 
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FIG. 6: Unfolded eigenvalue spacing distributions from mode 
family datasets TV = (137 values), N = 1 (214 values) and 
TV = 2 (178 values) calculated by CAS3D for our W7-X-like 
equilibrium, and the distribution for the combined spectrum, 
TV = 0, 1 and 2. 
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The W7-X variant equilibrium studied was generated 
with the VMEC [l^ code, which assumes the magnetic 
field to be integrable, so that all magnetic field lines lie on 
nested toroidal flux surfaces, which we label by s, the en- 
closed toroidal magnetic flux divided by the toroidal flux 
enclosed by the plasma boundary. The magnetic field is 
characterized on each flux surface by its winding number 
t(s). (In tokamaks its inverse, q = \/t, is more com- 
monly used.) As seen in Fig. 0] the rotational transform 
profile is nonmonotonic and has low shear (t ax i s = 1.1066, 
*min = 1-0491, *cd g o = 1-0754) so it is close to, but 
greater than, unity over the whole plasma. It is also seen 
from Fig. 0] equilibrium is interchange unstable, because 
the Mercier stability criterion is violated over the whole 
plasma. 

The CAS3D code expands the eigenfunctions in a fi- 
nite Fourier basis in the toroidal and poloidal angles, se- 
lected so as to include all \n\ < n max and all m such that 
n/m lies in a band including the range of *. The Fourier 
tableau is depicted graphically in Fig. [5] In this code, the 
radial dependence of the perturbation functions is treated 
by a hybrid Finite-Element approach, using a linear in- 



FIG. 7: Dyson-Mehta spectral rigidity as a function of subin- 
terval length L. Colour code: N=0 mode family: red; N=l: 
green; N=2: blue; combined data set: cyan. The Poisson- 
process limit is also indicated (black dashed). 



terpolation for the normal displacement and piecewise 
constant interpolations for the scalar components that 
describe the MHD displacement within the magnetic sur- 
faces. In the calculations discussed here, 301 radial grid 
points have been used. The kinetic energy was used as 
normalization, and therefore, the unstable eigenvalues A 
may be converted to a nondimensional growth rate 7 viz. 
7T A = #o(0)vTA|/£ o (0). Here, i? (0) is the major ra- 
dius and Bq(0) the equilibrium magnetic field measured 
on the magnetic axis. 

Because of the 5-fold symmetry of the equilibrium, any 
toroidal Fourier harmonic n in an cigenfunction is cou- 
pled to toroidal harmonics n ± 5. With the poloidal har- 
monics chosen to be positive, m > 0, there are just three 
uncoupled mode families TV = 0, 1, 2 (compare 2J). 
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FIG. 8: Unfolded eigenvalue spacing distributions for mode 
family N — 1, calculated using two different unfolding meth- 
ods. The results are seen to be consistent to within statistical 



We characterize the statistics of the ensembles of eigen- 
values within the three mode families using two standard 
measures from quantum chaos theory U Il6l Il7j , first 
renormalizing ("unfolding") the eigenvalues so their av- 
erage separation is unity. The first measure, shown in 
Fig. H3 is the probability distribution function P(x) for 
the eigenvalue separation x. The other, shown in Fig. 
is the Dyson-Mehta rigidity As(L), where L is the sub- 
range of unfolded eigenvalues used. 

As seen from Fig. when the statistics are analyzed 
within the three mode families the eigenvalue spacing dis- 
tribution function is closer to the Wigner conjecture form 
found for generic chaotic systems than to the Poisson 
distribution for separable systems, as might be expected 
from [l^. However, when the spectra from the three 
uncoupled mode familes are combined, there are enough 
accidental degeneracies that the spacing distribution be- 
comes close to Poissonian. 

To test the sensitivity to the precise method of un- 
folding chosen, we did the statistics using two different 
methods. They are the Gaussian unfolding an d a fit 
with exponentials. The results, shown in Fig. [Sj indicate 
little sensitivity to unfolding method. 



IV. CONCLUSION 

Although not presented here, when all unstable eigen- 
modes (i.e. all I, m, and n) are included, the eigenvalue 
spacing statistics for the ideal-MHD interchange eigen- 
value spectrum in a separable cylindrical approximation 
is close to that of generic separable wave equations 
despite our earlier finding 0,0 that the spectrum in the 
subspace of the most unstable radial eigenmode I = 
is nongeneric, as explained by the model presented in 
SecHU 

In this paper we have shown that a strongly three- 
dimensional stellarator equilibrium related to W7-X, 
the unstable interchange (Mercier) mode spectrum has, 
to within statistical uncertainties, similar statistics to 
generic quantum chaotic systems. That is, the over- 
whelming majority of eigenvalues are not "good quantum 
numbers" and can thus be expected to display sensitivity 
to small perturbations. This needs to be borne in mind 
when doing convergence studies using stability codes such 
as CAS3D. 

An interesting question for further work is whether 
other modes, such as drift waves, are quantum chaotic 
in stellarators, or if this is a pecularity of MHD modes. 
There is already evidence that kinetic effects make the 
semiclassical (WKB) dynamics closer to integrable [l9| . 

Another question is whether quantum chaos of Mercier 
modes occurs in machines with more field periods than 
the 5 of W7-X. Earlier work [2(J suggested that, in a 
10-field-period heliotron equilibrium related to the Large 
Helical Device (LHD), the spectrum is close to that of 
an equivalent axisymmetric torus, and thus not chaotic. 
However, a relatively few modes were studied and the 
spacing statistics were not calculated. 
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